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Abstract 


This final report summarizes the research carried out on the development of numerical 
approaches for simulating the plumes of electric propulsion devices. In Appendix I, we describe 
the development of a particle-based Monte Carlo approach for simulating the plumes and 
impingement of chemical propulsion systems. In Appendix II, we describe a plume model for an 
electric propulsion device (a Hall thruster). The research performed in this grant has led to 
significant advances in the capabilities for modeling these plume flows. 
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Overview 


The original premise of the proposed research was to develop a computer model for the 
plume of a particular spacecraft propulsion system (the arcjet) and to use the model to simulate 
synthesis of thin films. Early in the project, with the concurrence of the technical monitor, it was 
decided to concentrate on modeling the plume rather than the film synthesis processes. In 
addition, it was decided to focus not on arcjets but on the plumes of various chemical and electric 
propulsion devices. 

The plumes of spacecraft propulsion systems are of interest because of their potential 
deleterious effects on the host spacecraft. Before undertaking the present research, plume 
models of spacecraft propulsion systems were not well advanced. Simple, two-dimensional 
models existed for chemical rockets and semi-empirical models existed for Hall thrusters. The 
research performed under this grant significantly advanced the modeling capabilities for both 
these types of rocket engines. 

In terms of simulating the plumes of chemical rockets, the main activity involved the 
development of a three-dimensional, fully parallelized code employing the direct simulation 
Monte Carlo method (DSMC). The code was successfully verified against experimental data 
taken at DLR, Germany, for the three-dimensional impingement of a small rocket plume on a flat 
plate. The DSMC code was then applied to model impingement of a Primex hydrazine rocket 
plume on a typical spacecraft geometry. This work was published in the Journal of 
Thermophysics and Heat Transfer and the paper is included as Appendix I. 

A second investigation involved the development of a computer code for simulating the 
plumes of Hall thrusters which are an important form of electric propulsion. In this case, a 
combination of the DSMC and Particle In Cell (PIC) techniques was employed. The computer 
code was again parallelized. The model was applied to the SPT-100 Hall thruster and its 
physical accuracy assessed through direct comparison with several different experimental data 
sets. This work was published in the Journal of Spacecraft and Rockets and the paper is included 
as Appendix II. 

The research performed under this grant contributed significantly to two doctoral theses. 
One of the graduating students (Keith Kannenberg) now works at Lockheed Martin on spacecraft 
contamination, and the other (Doug VanGilder) works at Edwards AFB on the DSMC technique. 
Therefore, both the research and educational goals of the AASERT program have been 
successfully achieved. 
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Three-Dimensional Monte Carlo Simulations 
of Plume Impingement 


Keith C. Kannenberg* and Iain D. Boydf 
Cornell University, Ithaca, New York 14853-7501 


Three-dimensional plume impingement flows are simulated using the direct simulation Monte Carlo 
(DSMC) technique. TVo impingement problems are computed. The impact of a jet " f "'‘™ gen0 " c 
inclined flat plate is considered. Good agreement is found between surface quantities calculoted by DSMC 
and experimental data. A free molecular model of the plume is shown to provide reasonable estimates of 
impingement quantities. The plume of a hydrazine control thruster firing in a model satellite configuration 
is simulated. Surface quantities and net impingement effects are calculated. Free molecular analysis 
provides less accurate modeling of a multispecies flow with boundary-layer effects. 


Nomenclature 

A P = plume constant 

C = boundary-layer plume constant 

C P = specific heat at constant pressure 

p = impingement pressure 

q = heat flux 

R e = nozzle exit radius 

r* = throat radius 

T w = surface temperature 

T 0 , Po = stagnation temperature, pressure 

U, p = plume velocity, density 

a = local angle of incidence 

y = ratio of specific heats 

= boundary-layer thickness at exit 
r) t = net parallel efficiency 

0 lim = maximum inviscid turning angle 

0 O = boundary-layer streamline angle 

p * = throat density 

cr = accommodation coefficient 

r = shear stress 

Introduction 

S PACECRAFT in orbit require propulsion systems for var¬ 
ious functions such as attitude control and stationkeeping. 
Low-thrust rockets are often employed to meet these require¬ 
ments. The gas plume that is produced during a thruster firing 
may impinge on spacecraft surfaces. Impingement can have a 
number of undesirable effects such as the production of ad¬ 
ditional forces and torques, heat loads, and surface contami¬ 
nation. Each of these effects can reduce the overall lifespan of 
the spacecraft. Accurate modeling of plume impingement is an 
important factor in spacecraft design. 

Numerical simulations of thruster plumes can be used to 
characterize and predict impingement effects. At the extremely 
low densities seen in a plume expanding into vacuum, contin¬ 
uum fluid mechanics becomes invalid. The direct simulation 
Monte Carlo (DSMC) method simulates the gas at the micro¬ 
scopic level and can capture nonequilibrium effects occurring 
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mophysics and Heat Transfer Conference, Albuquerque, NM June 
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in a rarefied plume flow. DSMC has been successfully applied 
to nonequilibrium flows in propulsion, hypersonics, and ma¬ 
terials processing, 1 " 4 as well as to plume impingement flows. * 
The goal of this work is to validate the capability to simulate 
complex plume impingement flows using the DSMC tech¬ 
nique. To demonstrate this capability, two different three-di¬ 
mensional impingement flows are examined. The first flow 
considered is the impingement of a jet on an inclined flat plate. 
Experimental measurements of impingement quantities are 
used to verify the numerical results. The second flow is the 
impingement of a hydrazine thruster plume on a solar array in 
a model satellite configuration. Free molecular theory is used 
to model both flows analytically. 

Numerical Method 

The simulations presented in this study are performed using 
a three-dimensional version of MONACO, a numerically ef¬ 
ficient implementation of the DSMC method. MONACO is a 
general-purpose DSMC code designed for workstation archi¬ 
tectures. It is designed to be applicable to a wide range of flow 
problems without requiring problem-specific modifications. It 
can be run on a single machine or in parallel on an arbitrary 
number of processors. 

Free Molecular Analysis 

In the far field of a low-density plume expansion, the flow 
reaches very low densities and large Knudsen numbers. Under 
these highly rarefied conditions, the plume can be reasonably 
approximated as being collisionless. With this assumption, im¬ 
pingement quantities on a surface located within the plume can 
be calculated analytically. This free molecular analysis pro¬ 
vides a limiting case for surface quantities. 

Under the assumption of free molecular flow, surface quan¬ 
tities are a function only of the incoming ffeestream particles. 
Particles scattering off of the surface do not have any further 
interaction with the surface or the incoming gas. The net mass, 
momentum, and energy transferred to the surface is obtained 
by integrating over the incoming and outgoing particle distri¬ 
bution functions. 8 . 

If the incoming gas is assumed to have an equilibrium ve¬ 
locity distribution at the hypersonic limit, the pressure is given 
by 

p = ipU 2 { 2(2 - cr)(sin af + oVirKy- l)/y]V(T w /r 0 )sin a) 

Similar expressions are obtained for shear stress, heat transfer, 
and flux. 
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In the hypersonic limit, the fluid velocity is further approx¬ 
i mated b y the limiting result from isentropic theory, U hm = 

V2C P T» 

The density of the gas is modeled using the plume model 
of Simons. 9 This model describes the density as a function of 
the throat density and location in the plume expressed in polar 
coordinates measured from the center of the exit plane. The 
plume density is given by 

p(r, 0)/p* = A P (r*lr) 2 f(6) (2) 

m = {cos[(7r/2)(0/0 iim )]} 2/(Y " 1) , 0 < 0 O 

=/(0=ft,)<r w ‘*\ 00 <e<e liia (3) 

where 0 O is the angle between the plume axis and the stream* 
line separating the boundary layer from the isentropic core, 
and ft im is the maximum turning angle of a streamline at the 
exit for inviscid flow. The constant C is given by 

c = |A,V( y + My - 1)(^/25 £ ) <t - ,v<t+1) (4) 


These parameters are calculated analytically as a function of 
the stagnation conditions, gas properties, and source geome¬ 
ter 10 

try. 

The relations for plume density and velocity can be used 
with expressions for surface quantities such as Eq. (1) to obtain 
expressions for impingement quantities in terms of distance 
from the orifice, location in the plume, and the angle of inci¬ 
dence: 


—(?)W te) 

X j^2(2 - o)(sin a) 2 + a sin a j (5) 

/ \2 / y/(ri) / \ 

r=2op*A,m\j) ) cos « sin * 

q" = aV2RToPoA P f(0) [•j) -jj (6) 

I y T y y + l T w 1 . 

V y ~ 1 Lt " 1 2(7 - 1 ) T 0 J 


L Orifice - Plate Distance 
X Location on plate 
P Incidence Angle 
r Orifice exit radius 
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Fig. 1 Schematic of flat-plate impingement. 


method. The impingement surface is located a distance of 40 
mm from the orifice measured perpendicular to the surface 
(distance L in Fig. 1). The flow conditions and geometry yield 
a plume constant A P = 0.617 and a turning angle fti m = 130.5 
deg. 

Experimental Study 

The plate impingement configuration was investigated ex¬ 
perimentally at DLR, German Aerospace Research Center, 
Germany. Measurements of impingement pressure and shear 
stress were taken by Legge, 11 and measurements of heat flux 
were made by Dbring. 12 A variety of stagnation pressures, plate 
orientations, and separations were considered. 

Pressure and shear stress data were obtained using a pressure 
balance that directly measured the force on a floating element. 
Heat flux data were obtained by measuring the rate of change 
of surface temperature using thermocouples. Data were taken 
on the surface along a line that is coplanar with the plume 
axis. 

The experiments were performed in the high vacuum facility 
in Gottingen, Germany. The background pressure of the facil¬ 
ity was 0.045 Pa during pressure and shear measurements with 
the stagnation pressure considered in the present study (1000 
Pa). The background pressure was twice this value (0.090 Pa) 
during the heat transfer experiment. 

Pressure and shear-stress data were normalized to eliminate 
the effect of stagnation pressure and plate separation L. The 
following normalizations were used: 

P = (p/po)(L/r *) 2 (8) 


F=nU 


=& A„m 

m 




RTo 7-1 


(7) 


Inclined Flat Plate 

Impingement of a jet onto an inclined flat plate represents 
one of the simplest possible three-dimensional impingement 
problems. The relatively simple geometry allows it to be read¬ 
ily investigated both numerically and experimentally while still 
retaining three-dimensional effects. 

The problem under consideration is a freejet impacting on 
a flat plate. A plume of molecular nitrogen is generated by a 
sonic orifice. The orientation of the plate is varied relative to 
the axis of the plume. Figure 1 shows a schematic of the con¬ 
figuration. 

This study considers one particular set of orifice inlet 
conditions—unheated flow (T 0 = T w = 300 K) expanding from 
a stagnation pressure of 1000 Pa. The orifice itself is circular 
in cross section and 1 mm in radius, and the flow is assumed 
to be sonic at the exit. The Knudsen number based on orifice 
radius is 8 X 10‘ 3 at the orifice exit. These conditions ensure 
a rarefied plume flow suitable for calculation using the DSMC 


f=(r/p 0 )(L/r*) 2 (9) 

While this normalization does reduce the data to a significant 
degree, some dependence on the stagnation pressure is still 
observed. A reduction in stagnation pressure results in an in¬ 
crease in normalized pressure and shear stress as a result of 
the increasing rarefaction of the plume. For the present study, 
data obtained using a stagnation pressure of 1000 Pa are used 
whenever possible. 

Physical Modeling 

The inlet flow at the sonic orifice is modeled as a macro- 
scopically uniform stream of molecular nitrogen directed along 
the axis of the plume. The effects of a boundary layer at the 
exit are assumed to be small and are neglected. Inlet properties 
are calculated from isentropic theory assuming sonic condi¬ 
tions at the exit and stagnation conditions of 1000 Pa and 300 
K. This corresponds to a velocity of 323 m/s, a temperature 
of 252.2 K, and a number density of 1.53 X 10 23 m -3 . A 
background gas pressure is imposed to simulate the effects of 
the experimental tank pressure. 

The impingement surface is modeled assuming diffuse re¬ 
flection with full accommodation to a constant surface tern- 
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perature (a = 1.0). The surface temperature is equal to the 
stagnation temperature of the plume (300 K). 

Computational Modeling 

Fully unstructured grids using tetrahedral cells are used. 
Grids are generated using an advancing front package called 
FELISA. Cell sizes are approximately scaled to the local mean 
free path. Figure 2 shows an example of a grid for the case 
with the plume axis parallel to the surface (/3 = 0 deg). Only 
the boundary mesh is shown. The majority of cells are in the 
vicinity of the nozzle where the density is highest. Compres¬ 
sion at the surface is relatively small because of the low den¬ 
sities involved. As a result, the surface grid at the plate (not 
shown) consists of approximately uniform triangles for the 
cases considered. Symmetry planes are employed to reduce the 
size of the domain. Two planes are used in the normal im¬ 
pingement (P = 90 deg) case, and one each in the other cases. 

The simulation time step is varied across the domain, with 
each cell using a unique value. The local time step is scaled 
by the cell size that is approximately proportional to the in¬ 
verse of the local density. 

Simulations are load balanced for parallel execution using a 
simple scheme that groups cells according to the Z coordinate 
of their geometric center. The number of cells assigned to each 
processor is chosen to distribute particles evenly among pro¬ 
cessors. 

Topical parameters for the three-dimensional plume simu¬ 
lations are summarized next. Computational cost parameters 
for plate impingement: grid cells = 300,000; number of par¬ 
ticles = 3.5 million; transient steps = 22,500; sampling steps 
= 5,000; calculation time = 14 h; and steady state tj, = 95%. 
The calculations are performed using 16 nodes of an IBM 
SP-2. 

Comparison of Surface Quantities 

Impingement quantities are examined for several orienta¬ 
tions of the plate. Comparison with experimental data is used 
to verify the accuracy of the DSMC simulations of the plume 
and impingement. Comparison with predictions from the free 
molecular model provides an estimate of the accuracy of this 
theory for simple impingement flows as well as a qualitative 
basis for comparison with the simulation at the far ends of the 
surface where experimental data are not available. 

In the case of normal impingement of the plume (P = 90 
deg), the experimental data are taken along a radial line ex¬ 
tending out from the axis. The simulation data set includes 
data across the entire surface presented as a function of dis¬ 
tance from the axis. In the other cases, both the experimental 
and simulation data sets are taken along a line on the surface 
that is coplanar with the plume axis. 

Normal Impingement (fi m 90 Deg) 

Figure 3 compares normalized surface pressures for simu¬ 
lation, experiment, and free molecular theory. Good agreement 
is found between the data sets. The profiles show the same 



Fig. 2 Surface grid for p * 0-deg case. 


general shape: pressure falling from a maximum at the axis. 
The DSMC and theoretical profiles show an asymptote at the 
background pressure at large distances from the axis, where 
the plume has a minimal effect on the surface. 

The simulation underpredicts the experimental values by a 
significant amount (as much as 60%) at the far extent of the 
experimental data set. The smooth asymptotic shape seen in 
the DSMC and theoretical profiles is not as apparent in the 
measurements far from the axis. The slower decrease in pres¬ 
sure with distance from the axis may indicate difficulties in 
measuring pressures at the low densities away from the axis. 
It may also indicate a higher tank pressure in the facility than 
is reported. 

Normalized shear stress is considered in Fig. 4. Good agree¬ 
ment is observed between the simulation, experiment, and the¬ 
ory. The profiles each show the same shape with stress rising 
from zero at the plume axis to a maximum and then falling 
with decreasing gas pressure. DSMC stresses exceed the free 
molecular values far from the axis. This is likely a result of 
particles that strike the surface, undergo collisions with the 
incoming plume, and are scattered back to the surface at a 
high angle of attack. 

Surface heat flux is shown in Fig. 5. The simulation captures 
the trend of the experimental measurements and predicts the 
magnitude of the heat flux away from the axis. The simulation 
does overpredict the heat flux considerably near the axis (XIL 
< 0.75). 

Significant fluctuations appear in the numerical heat-flux re¬ 
sults, with the magnitude of the fluctuations increasing with 
distance from the axis. This is primarily a result of the sample 
size of particles striking the surface. TTie total number of par¬ 
ticles sampled per cell is on the order of 10 4 near the axis and 



Fig. 3 Comparison between DSMC and measured surface pres¬ 
sures at p m 90 deg. 



Fig. 4 Surface shear stress for p - 90 deg. Comparison between 
DSMC, experiment, and free molecular theory. 
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decreases with radial distance, up to a factor of 3 at the farthest 
radial extent. The statistical fluctuations in the DSMC simu¬ 
lation are magnified in the calculation of energy transfer that 
involves the second moment of the particle velocities. 

The theoretical curve has the same shape as the DSMC re¬ 
sult, but shows a larger magnitude of heat transfer. The dif¬ 
ference increases as the axis is approached. Gas densities are 
larger near the axis, and consequently, the mean free path and 
Knudsen number are smaller. The assumption of free molec¬ 
ular flow is weakest near the axis and collisional effects are 
more significant. This tends to lower the energy transfer below 
the free molecular value. The deviation is greater in terms of 
heat transfer than in pressure or shear stress because energy is 
a second moment of the particle distribution function. 

P - 45-Deg Case 

Figure 6 shows surface pressure profiles for the p = 45-deg 
case. Good agreement is agmh found between the simulation 
and experimental data sets as well as free molecular theory. 
The point of maximum pressure is just downstream of the 
point directly below the orifice, where the effects of density 
and angle of attack are high. The pressure drops off asymp¬ 
totically to the background value to both sides of this maxi¬ 
mum; this is more rapid behind the orifice. 

The experimental data are again higher than the DSMC and 
theoretical values at both ends of the profile. The magnitude 
of the pressures at these points are comparable to that seen at 
the end of the ft = 90-deg profile. This points toward uncer¬ 
tainty in pressure measurements at low densities or a higher 
back pressure in the facility. 



Fig. 5 Surface heat flux for p 85 90 deg. Comparison between 
DSMC, experiment, and free molecular theory. 



Fig. 6 Surface pressure for p ■ 45 deg. Comparison between 
DSMC, experiment, and free molecular theory. 



Fig. 7 Surface shear stress for p - 45 deg. Comparison between 
DSMC, experiment, and free molecular theory. 



Fig. 8 Surface heat flux for p « 45 deg. Comparison between 
DSMC, experiment, and free molecular theory. 


Shear-stress profiles are reported in Fig. 7. The DSMC 
stresses show excellent agreement with the experimental data. 
The features of the experimental profile are captured by the 
simulation, particularly the rapid decrease that occurs directly 
below the orifice. 

While the qualitative structure of the free molecular profile 
agrees with those of the simulation and experiment, the theo¬ 
retical profile differs in terms of magnitude. In front of the 
orifice (X/L > 0), the free molecular stresses overpredict the 
simulation and measured values. Behind the orifice (X/L < 0), 
the free molecular stresses underpredict the actual stresses. 
This is likely a result of the presence of backscattered particles 
that strike the surface behind the orifice. 

Heat transfer profiles are reported in Fig. 8. Although sig¬ 
nificant fluctuations are again seen in the simulation and ex¬ 
perimental results, there is generally good agreement between 
the data sets. Free molecular theory overpredicts heat transfer 
across most of the range considered because of collisional ef¬ 
fects. This effect is countered to some extent behind the orifice 
by backscattered particles that tend to increase the transport of 
energy to the surface. 

Parallel Impingement (p - 0 Deg) 

Before presenting comparisons with experimental data, it is 
important to note that the DLR study by Legge" does not 
report pressure or shear-stress data for the parallel impinge¬ 
ment case with a stagnation pressure of 1000 Pa. The lowest 
stagnation pressure for which J3 = 0-deg pressure and shear 
data are reported is 4000 Pa. The study by Doring 12 does pre¬ 
sent heat flux data for a stagnation pressure of 1000 Pa. The 
DSMC simulations use a stagnation pressure of 1000 Pa. Al¬ 
though the data are normalized by the stagnation pressure. 
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Legge reports some effect of flow rarefaction on the normal¬ 
ized surface quantities, particularly shear stress. 

The difference in stagnation pressures complicates the mod¬ 
eling of tank pressure effects. The ratio of background to stag¬ 
nation pressure is not a constant in the experimental facility, 
generally falling with increasing stagnation pressure. Conse¬ 
quently, the parallel impingement case has been simulated 
without backpressure and the theoretical pressure profile is for 
expansion into a vacuum. 

Surface pressures for the parallel impingement case are 
shown in Fig. 9. The profiles show good qualitative agreement. 
Without a backpressure, the simulation and theoretical profiles 
drop toward zero at the ends of the plate. No asymptote is seen 



Fig. 12 Schematic of model satellite configuration. 



Fig. 9 Surface pressure for P - 0 deg. Comparison between 
DSMC, experiment, and free molecular theory. 



Fig. 10 Surface shear stress for 0 * 0 deg. Comparison between 
DSMC, experiment, and free molecular theory. 



Fig. 11 Surface heat flux for 0 * 0 deg. Comparison between 
DSMC, experiment, and free molecular theory. 


in the experimental data because the normalized facility pres¬ 
sure (not shown) is significantly lower with the higher stag¬ 
nation pressure. The lower experimental pressures are a result 
of a higher stagnation pressure and lower Knudsen number at 
the surface. 

Shear stress data are shown in Fig. 10. The profiles again 
show qualitative agreement. The experimental values are lower 
because of higher densities. The effect on shear is stronger 
than that on pressure. This was previously reported by Legge. 
The effect of rarefaction can be seen to a much smaller degree 
in a comparison between simulation and free molecular results. 
The DSMC values are slightly lower because of collisional 
effects. 

Heat transfer results are presented in Fig. 11. Experimental 
heat-flux data are based on a stagnation pressure of 1000 Pa, 
equal to that used in the DSMC simulation. At the low den¬ 
sities involved in the parallel impingement case, statistical fluc¬ 
tuations are quite large in both the experimental and numerical 
results. The simulation results overpredict the data to some 
degree, although the trend is captured quite well. The free 
molecular results do not significantly overpredict the heat 
transfer as was seen in the other cases. This is likely a result 
of the low densities found at the surface for this case. 

Model Satellite Configuration 

A representative satellite geometry consisting of a spacecraft 
bus and solar arrays is simulated as a model impingement con¬ 
figuration. A plume is generated by a hydrazine control thruster 
mounted at the comer of one side of the spacecraft bus. Fig. 
12 shows a schematic of the problem. The spacecraft bus is a 
1.5-m cube of which half is simulated. The array section is 
3.25 m long, 2.5 m wide, and 0.12 m thick, and is deployed 
0.75 m from the spacecraft. Expansion of the plume and its 
impingement on the solar array panel are modeled. Two ori¬ 
entations of the array with respect to the bus are considered. 

The thruster is modeled after an operating hydrazine control 
thruster, the MR-103-series 0.2-lbf REA, manufactured by Pri- 
mex Aerospace Company. This thruster was originally devel¬ 
oped for the Voyager 1 and 2 spacecraft and is also in use on 
such programs as GPS, Iridium, Cassini, and various com¬ 
munications satellites. This study will consider an operating 
mode at 60% of the rated thrust. Relevant data for the thruster$ 
are shown next: Operating conditions and parameters for hy¬ 
drazine thruster: nominal thrust — 0.55 N, expansion ratio — 
100:1, half angle = 15 deg, exit radius = 2.921 mm, chamber 
pressure = 1.25 MFa, stagnation temperature = 1167 K, and 
flow rate = 2.62 X 10“ 4 kg/s. 

Physical Modeling 

The gas plume is a three-species mixture of molecular ni¬ 
trogen, hydrogen, and ammonia. Simulations begin at the exit 
plane of thruster using flow properties based on a numerical 
calculation of flow in the thruster. Exit plane number densities 

^Morris, J., private communication, Primex Aerospace Co., Red¬ 
man, WA, 1997. 
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taken from this earlier calculation are scaled down by a factor 
of 1.52 to maintain consistency with the experimentally mea¬ 
sured mass flow of the device. At the centerline the inlet flow 
has the following properties: 2200 m/s velocity, 500 K tem¬ 
perature, and 3.2 X 10 23 m~ 3 total density. 

Because of the relatively low temperature of the plume, it 
can be considered chemically frozen. Vibrational energy 
modes are likewise frozen. All spacecraft surfaces are modeled 
as diffuse reflectors with full accommodation. Surfaces are 
maintained at 273 K. Surface chemistry and adsorption are not 
considered in these calculations. 

Computational Modeling 

Simulations are performed in two parts. First, an axisym- 
metric calculation of the near-field plume in the vicinity of the 
thruster is performed. Data from this calculation are then used 
as input for the three-dimensional calculation of the plume far 
field and impingement. Separation of the calculation in this 
manner greatly reduces the cost of the overall calculation. 

Near Field 

The near field expansion of the plume is considered to be 
an axisymmetric problem that is independent of the surround¬ 
ing geometry. The flow domain is hemispherical and extends 
more than 10 exit radii out from the center of the thruster exit. 
The grid is an unstructured mesh with cell sizes and time steps 
scaled according to the local mean free path. 

Inflow conditions for the far-field calculation are taken at 
the breakdown surface. This is the locus of points for which 
the breakdown parameter 13 equals the threshold value P = 0.02. 
To more readily use this as a geometric interface between the 
near- and far-field simulations, the surface is simplified to a 
capped cylinder with uniform cross section. Density in the near 



Fig. 13 Contours of number density for the axisymmetric near¬ 
field simulation. The interface surface is indicated. 



Fig. 14 Boundaries of the computational domain. Shown for un¬ 
rotated array case, * 0 deg. 


field is shown in Fig. 13 along with the location of the inter¬ 
face surface. 

Far Field and Impingement 

The far-field simulation considers the whole area surround¬ 
ing the spacecraft, including the solar array. The computational 
domain is an orthogonal parallelepiped that extends beyond 
the spacecraft surface by at least 0.25 m in each coordinate 
direction. Only one-half of the spacecraft is simulated. The 
boundaries of the computational domain are shown in Fig. 14. 
The solar array is shown in the nominal, unrotated configu¬ 
ration (<F = 0 deg). 

An unstructured mesh consisting of tetrahedral cells is em¬ 
ployed. The mesh is generated using the grid-generation pack¬ 
age FELISA. Cell sizes are approximately scaled according to 
local density. Variable time steps are employed with cell size 
being used as a scaling factor. 

To improve resolution in the vicinity of the interface surface 
where flow density is highest, a simple scheme of particle 
weight scaling is employed. Cells in the immediate vicinity of 
the inflow (interface) surface are assigned a low particle 
weight, 10% of the reference value for the simulation. The 
band of cells immediately outside this surface are assigned a 
weight of 40% of the reference value. The remaining cells are 
assigned the reference weight. Figure 15 shows this weight 
scheme on a plane cut taken along the axis of the cylindrical 
inflow surface. This scaling roughly follows the density vari¬ 
ation in close proximity to the inflow. Because there is only a 
small number of particles moving toward the inflow in this 
highly supersonic flow, there is a minimal amount of cloning 
involved. The weight scheme allows reasonable resolution of 
the high-density region without affecting the accuracy of the 
calculation. 

The computational parameters for the near- and far-field 
plume calculations are summarized in Table 1. 

Results 

Impingement simulations were performed for two orienta¬ 
tions of the solar array, = 0 deg and <I> = +30 deg. Of 


'Rtble 1 ^Computational cost parameters for model 
satellite configuration 


Parameter 

Near field 

Far field 

Grid size 

7600 cells 

290,000 cells 

Number of particles 

425,000 

2,900,000 

Transient steps 

12,500 

18,000 

Sampling steps 

10,000 

10,000 

Calculation time 

5 h 

11 h 

Number of processors 

4 

16 

Parallel efficiency 

97% 

91% 



Fig. 15 Particle weight scaling in the vicinity of the inflow sur¬ 
face. 
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primary importance are the mass, momentum, and energy im¬ 
parted to the array by the plume. Figures 16 and 17 show 
distributions of surface pressure for the 0- and +30-deg cases. 
The view is in the plane of the array. The contour plots show 
maximum values at the near end of the array toward the +T 
side, closest to the thruster. In the <b = 0-deg case, the maxi¬ 
mum pressure is downstream of the leading edge as a result 
of the influence of the nozzle boundary layer. As expected, 
impingement pressures are larger in the 4> = +30-deg case, 
where the array is tilted toward the thruster. 

Ammonia molecules are the major contamination concern 
with hydrazine thrusters. Figures 18 and 19 show distributions 
of ammonia flux on the array. The qualitative behavior of this 
property is very similar to that observed with pressure. 



Fig. 16 Contours of impingement pressure (Pa) at the array sur¬ 
face for the 0-deg array case. 



Fig. 17 Contours of impingement pressure (Pa) at the array sur¬ 
face for the 30-deg array case. 



Fig. 18 Contours of NH 3 flux (no./m 2 s) at the array surface for 
the 0-deg array case. 


Data extracted along a line on the array surface are used to 
quantitatively examine the surface quantities. Data are taken 
on the line that intersects the array top surface and the plane 
parallel to the array axis that contains the plume axis. 

The operating conditions for the hydrazine thruster can be 
used to calculate the plume parameters A P > 6 0t and C. 
There is some uncertainty involved in this calculation when a 
multispecies gas is considered. To be most useful from an en¬ 
gineering standpoint, the free molecular results should be 
based on stagnation conditions and the geometry of the nozzle. 
The expressions for surface quantities [Eqs. (5-7)] depend on 
gas properties R and y, which are difficult to determine a priori 
for a rarefied gas in which individual species will have varying 
mole fractions and velocities. In this study, the gas has been 



Fig. 19 Contours of NH 3 flux (no./m 2 s) at the array surface for 
the 30-deg array case. 



Distance Along Array (m) 

Fig. 20 Comparison of simulation and free molecular pressures 
on array surface for 0-deg array case. 



Distance Along Array (m) 

Fig. 21 Comparison of simulation and free molecular pressures 
on array surface for 30-deg array case. 
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assumed to have a constant molecular weight and ratio of spe¬ 
cific heats using values reported for the exit plane (M =13 
kg/mol, y = 1.4). These parameters yield a plume constant, A P 
= 3.62, boundary-layer constant, C = 12.6, and limiting angle, 
0 O = 19.2 deg and 4m = 58.7 deg. 

Figure 20 shows impingement pressure profiles for the 0- 
deg case. Reasonable agreement is found on the far end of the 
array, but at the near end the free molecular result falls well 
below the simulation. The noticeable drop in the free molec¬ 
ular profile on the upstream side of the array is a result of 
boundary-layer effects. The drop in pressure can be seen in 
the simulation as well, but it is not as dramatic. The sharp 
decline indicates that the density model falls too rapidly in the 
boundary-layer regime. 



Fig. 22 Comparison of simulation and free molecular heat flux 
on array surface for 0-deg array case. 



Fig. 23 Comparison of simulation and free molecular heat flux 
on array surface for 30-deg array case. 



Fig. 24 Comparison of simulation and free molecular NH 3 flux 
on array surface for 0-deg array case. 


The 30-deg case is shown in Fig. 21. In this case, the array 
is rotated toward the thruster, and as a result the profile being 
considered is significantly closer to the thruster. The linear 
profile is within the core of the plume, and consequently, there 
is no drop in pressure on the upstream side of the plate as a 
result of the boundary layer. Both the simulation and analytical 
profiles primarily show an inverse squared drop in pressure 
with distance, which is consistent with the density model. The 
DSMC results predict a 70% higher pressure across the length 
of the array. This may indicate a low estimate of the plume 
constant A P . 



Fig. 25 Comparison of simulation and free molecular NH 3 flux 
on array surface for 30-deg array case. 



Fig. 26 Comparison of NH 3 densities from stand-alone plume 
simulation and free molecular model. Data at location of array 
surface for 0-deg case. 



Fig, 27 Comparison of NH 3 velocities from stand-alone plume 
simulation and free molecular model. Data at location of array 
surface for 0-deg case. 
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Figures 22 and 23 show heat transfer profiles. In both cases 
the behavior is very similar to that of surface pressure. The 
free molecular profile for 0 deg shows an extreme drop on the 
upstream end while otherwise matching the shape of the 
DSMC profile. In the 30-deg case the DSMC and theoretical 
profiles are qualitatively similar. In both cases the free molec¬ 
ular values tend to exceed the simulation heat transfer. Over- 
prediction of heat flux is consistent with the results for flat- 
plate impingement. 

The number density flux of ammonia striking the array is 
shown in Figs. 24 and 25. Two free molecular profiles are 
shown for each case. The first (labeled A) assumes a constant 
species composition in the plume model and a uniform limiting 
velocity based on the average molecular weight. The second 
(B) calculates the ammonia flux as if the plume was composed 
entirely of ammonia with the correct number density. The two 
methods produce essentially the same fluxes. 

Qualitatively, thp flux profiles show the same behavior for 
both pressure and heat transfer. Both angle cases show signif¬ 
icantly higher DSMC fluxes. This indicates the difficulty of 
applying the free molecular model to a multispecies, rarefied 
gas. 

The problems with the free molecular model for flux are 
illustrated by a simulation of the thruster plume flowfield with¬ 
out any impingement surfaces. This plume flow should give a 
reasonable estimate of what the impingement surfaces see, be¬ 
cause the boundary layer caused by the impingement is small. 
Figure 26 plots ammonia density for the freely expanding 
plume at the same spatial location as the linear profiles shown 
for the O-deg case in Fig. 24. The DSMC results are compared 
with density predicted by the Simons plume model (assuming 
an ammonia plume, method B discussed earlier). 

Figure 27 compares simulated ammonia velocities and the 
li miting velocity employed in the free molecular analysis 
(\Z2C P T 0 ). Downstream of where the leading edge of the solar 
array would be located, within the isentropic core of the plume, 
DSMC densities are on average two times higher than the Si¬ 
mons model predicts. DSMC velocities are also 5-8% larger 


Table 2 Integrated impingement effects for * 0 deg 


Impingement property 

DSMC 

Free molecular 

Net force, X component, N 

0.100 

0.118 

Net force, Y component, N 

-0.0107 

-0.00588 

Net force, Z component, N 

-0.0814 

-0.0647 

Net centerline torque, Nm 

0.0405 

0.0369 

Net heat transfer, W 

101 

112 

Net incidence of NH 3 , s“* 

1.44 X 10 21 

470 X 10“ 


Table 3 Integrated impingement effects for 

■ +30 deg 

Impingement property 

DSMC 

Free molecular 

Net force, X component, N 

0.0922 

0.148 

Net force, Y component, N 

0.0228 

0.0226 

Net force, Z component, N 

-0.0645 

-0.0563 

Net centerline torque, Nm 

0.0519 

0.0509 

Net heat transfer, W 

88.3 

137 

Net incidence of NH 3 , s" 1 

1.43 X 10 21 

6.15 X 10 20 


than the analytical model. Because the number density flux in 
free molecular flow is the product of density and velocity [Eq. 
(7)], these factors combine to produce the significantly larger 
fluxes seen in the DSMC results shown in Fig. 24. 

In the boundary-layer region, the relationship between 
plume properties and impingement flux is less clear. However, 
near the upstream edge, the simulation again produces signif¬ 
icantly higher ammonia densities, which again would lead to 
larger flux. The more complex density profile seen in the 
DSMC result, likely caused by differing mole fractions and 
velocity slip between the species, again indicates the limita¬ 
tions of the analytical plume model for complex gas flows. 

The total forces and torques imparted to the array are im¬ 
portant for spacecraft design. The total amount of a contami¬ 
nant such as ammonia striking the surface is also important. 
Tables 2 and 3 summarize these integrated quantities for the 
two simulations. Integrated quantities calculated from free mo¬ 
lecular theoiy are also included. 

It is worthy to note that in both cases the net force exerted 
on the array by the plume is on the order of 20% of the nom¬ 
inal 0.55 N thrust of the thruster. This interaction force is pri¬ 
marily in the +X direction, which is opposite to the thrust 
vector and thus acts to reduce the effective thrust of the device. 
As a result, nearly 20% of the propulsion energy is wasted. 

Little net difference is seen between the two array orienta¬ 
tions. Although part of the array is closer to the thruster in the 
30-deg case, the average angle of attack over the array is 
higher, which reduces the impingement quantities. The rotation 
causes a portion of the array to be significantly farther from 
the thruster, reducing the effective area over which force is 
exerted. 

Comparison between the simulation and free molecular val¬ 
ues follows the same pattern shown in the linear profiles. The 
net forces are comparable, which is consistent with similar 
pressures. Free molecular flow predicts a higher energy trans¬ 
fer and a significantly lower net flux of ammonia, both con¬ 
sistent with the previous comparisons. 

In all of the previous analyses and simulations, surfaces 
were assumed to be fully accommodating. A more accurate 
representation of a real surface would assume partial accom¬ 
modation. An accommodation coefficient a of 0.8 is typical 
for a metallic surface. Table 4 compares integrated impinge¬ 
ment effects for full and partial accommodation. As expected, 
net heat transfer to the array scales with accommodation co¬ 
efficient. The net force in the x direction, parallel to the plume 
axis, is a result of surface shear, and consequently, also scales 
with or. The other two force components and torque are essen¬ 
tially unaffected by a change in accommodation coefficient. 
These quantities are primarily caused by surface pressure that 
scales nonlinearly with cr. 

Examination of the comparisons in Table 4 and the expres¬ 
sion for impingement pressure [Eq. (5)] indicate that the sur¬ 
face model will have some effect on the net impingement ef¬ 
fects. Heat transfer will be overpredicted by a fully diffuse 
model. Momentum transfer will be overpredicted if it is pri¬ 
marily a result of shear forces. If more direct impingement is 
involved, such as might occur in a docking maneuver, the ef¬ 
fect is less clear and will depend on the surface and plume 
temperatures and the geometry through angle of attack. These 


Table 4 Effect of accommodation coefficient on free molecular integrated impingement effects 


Impingement property 

$ = 0 deg, 
cr= 1.0 

4> - 0 deg, 
cr — 0.8 

^ = +30 deg, 
a = 1.0 

= +30 deg, 
o-= 0.8 

Net force, X component, N 

0.118 

0.0945 

0.148 

0.118 

Net force, Y component, N 

-0.00588 

-0.00471 

0.0226 

0.0227 

Net force, Z component, N 

-0.0647 

-0.0656 

MMiBinmwr.ic^Bli 

1 1 1 m mm 

Net centerline torque, Nm 

0.0369 

0.0377 

0.0509 

0.0490 

Net heat transfer, W 

112 

89.5 

137 

109 
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observations should hold for both an analytical model and 
DSMC simulations. 

Conclusions 

The DSMC method provides reasonably accurate modeling 
of impingement flows. Generally good agreement with ex¬ 
perimental data for the plate problem indicates that the plume 
and impingement effects are modeled correctly, although there 
is some tendency to overpredict surface heat flux. The calcu¬ 
lation of a satellite configuration demonstrates the ability to 
simulate real engineering configurations. Calculations are ex¬ 
pensive in three dimensions, but can be performed efficiently 
in parallel. Careful use of variable scaling can reduce the over¬ 
all cost significantly. 

Free molecular theory provides a reasonable estimate of sur¬ 
face quantities at high Knudsen numbers. The analysis tends 
to overpredict the values, however, with heat transfer and shear 
stress being more sensitive to the degree of rarefaction than 
pressure. Boundary-layer and multispecies effects are not prop¬ 
erly handled by the simple plume model. 

Problems that can be broken down into two distinct parts 
such as the near and far field of a plume can be more effi¬ 
ciently computed using two separate simulations. If the first 
portion can be computed assuming axial symmetry, the overall 
cost can be reduced significantly. This type of hybrid approach 
is particularly appropriate for engineering problems, where 
simulations of a number of far-field problems may be begun 
using a single near-field result. 
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Particle Simulations of a Hall Thruster Plume 


Douglas B. VanGilder,* Iain D. Boyd* and Michael Keidar* 
Cornell University, Ithaca, New York 14853 


A numerical code that combines the direct simulation Monte Carlo method with the particle-in-cell method 
has been developed to examine the plumes of Hall thrusters. The present investigation includes a study of the 
sensitivity of the computed plume to various ion conditions at the thruster exit and considers models for computing 
the electric field based on the electron momentum equation. Specifically, two electrostatic models are compared: 
one assumes isothermal electrons, whereas the other uses a variable electron temperature model. Computations 
are compared with experimental measurements of current density, ion velocity, ion density, electron density, heat 
flux, and plasma potential in the plume of a stationary plasma thruster. The varying electron temperature is found 
to affect the very near field significantly. Simulations using this model agree better with near-field current density 
measurements. This model also leads to better agreement with electron number density measurements in the 
far field. The agreement of plasma potential depends on parameters used in the simulations. To better represent 
the flow inside the experimental facilities, the full chamber geometry is simulated assuming symmetry about the 
thruster centerline. These simulations lead to excellent agreement with current density measurements to 180 deg. 
The sensitivity study indicates that the far-fieid current density is relatively easy to reproduce. 


Nomenclature 

E = electric field, V/m 
e = charge of an electron, C 
k = Boltzmann constant, J/K 
nj = number density of species j, m~ 3 
n n f = reference electron number density, m -3 ; defined where 
potential is 0 
p = pressure, N/m 2 

T e = electron temperature, K unless specified as eV 
Wj - particle weight of species j 

Xe = xenon 

y = ratio of specific heats 

<f> = plasma potential, V 

Subscripts 

e = electrons 

/ = fast (high energy) 

/ = ions 

n = neutrals 

r = radial direction 

s = slow (low energy) 

x = axial direction 

Superscript 

+ = ion, either singly or doubly charged 

Introduction 

ITH increased commercial applications of satellites, there is 
a need for onboard thrusters with relatively long lifetimes. 
As with other electric propulsion devices, Hall thrusters offer a high 
specific impulse that is well suited for satellite station-keeping, repo¬ 
sitioning, and orbit transfer. The stationary plasma thruster (SPT) 
variety has been used for many years. Xenon gas is currently the 
propellent of choice for such devices because it is an inert gas with a 
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relatively low ionization potential. There is concern, however, about 
contamination of the satellite surface caused by the plasma in the 
plume. The charge-exchange plasma, created by collisions between 
ions and nonionized propellent in which electrons are transferred, is 
of particular concern. The charge-exchange ions have much lower 
velocities than the propellent ions; therefore, they are more influ¬ 
enced by the self-consistent electric fields. These fields may cause 
them to interact with spacecraft surfaces. The charged particles may 
impact solar arrays or interfere with transmissions from the satellite. 
Therefore, it is important to understand in detail the dynamics of 
the plume. 

Computational modeling allows the dynamics of the plume and 
interaction with its environment to be examined without any influ¬ 
ences of experimental facilities. The ability to simulate the plumes 
of these devices allows a wider variety of operating conditions to be 
tested than would be feasible in a laboratory. Many test chambers 
can be maintained at fractions of millipascals. However, the oper¬ 
ating conditions of these devices give relatively low densities, and 
thus interactions between the plasma and the background gas are 
to be expected. Also, the walls of the chamber may influence the 
plasma. However, computer codes must rely on experimental data 
for validation. The goal of the present study is to assess the ability 
to simulate these plumes accurately. 

To accomplish this task, a computer code that combines the direct 
simulation Monte Carlo 1 (DSMC) and the particle-in-cell 2 (PIC) 
techniques is being developed to understand in detail the plasma 
behavior of the plumes of Hall thrusters. The PIC method deter¬ 
mines the trajectories of charged particles as predicted by imposed 
and self-consistent electric fields. The DSMC method is used to sim¬ 
ulate the collisional effects in the flowfield. Both charge-exchange 
and momentum-transfer collisions are modeled. Ions and neutral 
atoms from the thruster and background atoms are simulated. The 
code has previously been verified for an ion thruster 3 and is being 
applied to Hall thrusters in the current study. The plumes of these 
two types of thrusters are similar. The main differences in modeling 
are in the geometry and the ratio of neutral atoms to ions. Numerical 
improvements to this earlier version of the code are discussed in the 
next section. 

Because electric propulsion devices such as ion thrusters, arc jets, 
and Hall thrusters inherently involve charged propellent, the PIC 
technique is well suited to simulate their plumes^ Roy 4 developed 
a code for modeling ion thrusters. Wang et al. 5 have also modeled 
ion thrusters. Previous work by Oh and Hastings 6,7 modeled SPT 
plumes using a PIC-DSMC code. In contrast to the work presented 
in Refs. 4 and 5, the neutrals are also tracked in the present study, 
and charge-exchange ions are generated directly from collisions 
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between ions and neutrals. The underlying code (MONACO ) is a 
parallel DSMC code to which PIC routines are added. Extensive 
experimental data have been obtained for the SPT- 100 thruster con¬ 
sidered in the present study. 9 -' 4 Measurements of current density 
and ion velocity were reported in Refs. 10 and 11. Recently, mea- 
surements of ion velocity, ion density, heat flux, neutral flux, and 
plasma potential were obtained. 1213 (The reader is also referred to 
a private communication by C. Marrese and A, Gallimore, Uni¬ 
versity of Michigan, May 1999.) The present study relies on these 
new data to perform the most comprehensive comparisons between 
experimental and computational work thus far. 

The following subsections describe the assumptions and the nu¬ 
merical approach to the modeling, including a brief explanation 
of the numerical techniques used. The various simulations are de¬ 
scribed, followed by a discussion of the results and some conclu¬ 
sions. 

Physical Modeling 

The plume is comprised of propellent ions, charge-exchange ions, 
neutral atoms from the thruster, electrons, and the background gas 
of the experimental facility. The interactions of these species as well 
as the influence of the electric and magnetic fields are the important 
modeling issues. The approach to this modeling is outlined m the 

following paragraphs. 4 . 

Although utilization efficiencies are high for Hall thrusters, the 
nonionized propellent accounts for about 50% of the density at the 
exit. The neutral density is substantially higher when the cathode 
flow is included. The neutrals may undergo collisions with the ions. 
Many of these collisions are charge-exchange reactions, where an 
electron is transferred from a neutral at a relatively low velocity to 
a high-energy ion. The probability of such a collision is a function 
of relative velocity and is modeled using the cross section given by 
Rapp and Francis. 15 This process creates the charge-exchange ions 
already mentioned. At the thruster exit the nonionized propellent is 
modeled at sonic conditions based on a temperature of 1000 K. These 
assumptions yield a neutral velocity of 325 m/s, which agrees well 
with experiments with a similar Hall thruster. 16 The sonic assump¬ 
tion has been shown to lead to good agreement with experiments of 

neutral xenon flow. 17 . 

The background gas of laboratory facilities may interact with the 
propellent during experiments. Although the experiments consid¬ 
ered here were conducteu with a back pressure of 6 mPa, at room 
temperature this gives a density on the order of the thruster exit 
neutral density. The background gas is assumed to be composed en¬ 
tirely of xenon neutrals at 300 K and to be uniform. Collisions with 
the species that originated from the thruster are simulated, but the 
background particles are not tracked by the code. Instead, tempo¬ 
rary particles are created in each cell with velocities chosen from a 
Maxwellian distribution and are paired for possible collisions with 
the species from the thruster. These may be either momentum trans¬ 
fer or charge-exchange reactions. At distances of more than a few 
centimeters from the exit plane, the background neufrals become 
the dominant source of charge-exchange ion production. As with 
the ion thruster plume considered in a previous study, the changes 
to the neutral distribution function would be due to the fast atoms 
whose density would be at least three orders of magnitude lower 
than the bulk density. Thus, it is reasonable to assume that the back¬ 
ground gas is unchanged with collisions, and it is unnecessary to 
use computational time and memory on them. 

As with the previous study, 3 the ions are tracked, but the electrons 
are not. Instead, the electron behavior is given by a balance of the 
electrostatic force and the pressure gradient. Mathematically, this is 

n e eE = —Vp 

The balance is necessary to prevent the electrons from leaving a re¬ 
gion en masse, which would leave behind a large charge imbalance. 
This equation comes from the electron momentum equation and ne¬ 
glects the collision term and the magnetic field. For the densities and 
temperatures of the plasma typical of the plumes of Hall thrusters, 
the ion-electron collision term is negligible. The ratio of the col¬ 
lision frequency to the plasma frequency is much less than 1. The 
magnetic field term is generally neglected. 3 ’ 4,6 This appears to be 


a reasonable estimate for the far field. Further assuming isothermal 
conditions yields the familiar Boltzmann relation 

n e = n ref exp(e^/ kT e ) 

However, the near field in the plumes of Hall thrusters has variations 
in the electron temperature. The effects of including this variation 
are considered in the present study. For simplicity, the potential 
is not calculated in the simulations. Instead, the electric field in the 
axial and radial direction are calculated directly. Assuming adiabatic 
conditions, this becomes 

kT t db(n e ) k d T, 

Exr ~~ e d(jc, r) ed(x,r) 

An electron temperature is imposed in the domain based on mea¬ 
surements taken by Kim et al. 15 The electrons are assumed to act as 
an expanding fluid at isentropic conditions. Their density decays as 
r” 2 and the temperature scales with n Y t . Constants are chosen to 
closely match the measurements of Kim. A comparison in Fig. 1 in¬ 
dicates the good agreement obtained with this approach. A far-field 
temperature is set for the simulations as well. Simulations are run 
using both models to determine the effects on the flowfield. 

The difference in magnitudes of the ion and electron velocities 
makes it difficult to also track the electrons in the simulation. By 
assuming that the plasma is quasineutral (it/ *n e ), the ion density 
can be used to find the potential. Both singly and doubly charged 
propellent ions as well as the charge-exchange ions are included m 
the simulations. Various ion inlet conditions are used for these sim¬ 
ulations to determine whether the comparisons with experimental 
data are sensitive to the ion profile at the thruster exit. 

Figure 2 illustrates the computational domain used for most of 
the simulations. The domain is assumed to be axisymmetric about 
the thruster centerline. The main three-dimensional effects would 
come from the neutralizing cathode. The flow through the cathode is 



Fig. 1 Comparison of analytical electron temperature to measure¬ 
ments by Kim. 14 
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Fig. 2 Computational domain. 
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assumed to exit from the thruster as neutral xenon atoms. The radial 
electric field is set to zero on the centerline to satisfy symmetry, 
but axial variations in potential are permitted. The boundary at the 
maximum radial position of the domain has the same conditions. 
The two axial boundaries permit only radial variations in potential. 
For an adequately large domain the variations at the boundaries 
would be negligible. Particles that reach boundaries other than the 
symmetry line leave the simulation. The geometry below the exit 
plane is simplified. A wall at ground potential is assumed to extend 
from this exit to the centerline. When an ion strikes a wall, it is 
neutralized. Thus the assumption is made that an electron strikes the 
wall as well to maintain quasi-neutrality. The thruster wall boundary 
condition has little effect on the flow values compared in the results 
section. 

Larger simulations include the full geometry of the experimen¬ 
tal facility used at the University of Michigan as the computational 
domain to better represent the experimental conditions 13 in the sim¬ 
ulations. This domain has a length of 9 m and a radius of 3 m. The 
walls are assumed to be at ground potential. Particles that reach the 
walls are removed because the background routine already main¬ 
tains the facility background pressure. Again the flow is assumed to 
be symmetric about the thruster centerline. 

Computational Modeling 

Because both the DSMC and PIC numerical methods rely on par¬ 
ticles that represent real molecules, it is desirable to have enough 
particles to represent the velocity distribution adequately. The fol¬ 
lowing techniques permit larger simulations than would otherwise 
be feasible. 

The underlying DSMC code 8 is parallelized to run on the IBM 
SP/SP2 architecture. The computational cells are divided among 
the processors such that each processor maintains its own collection 
of neighboring cells. Parallel PIC codes also partition by blocks of 
cells. 18 However, in the present implementation separate grids are 
maintained by the two algorithms. The data structures are defined 
such that the particles are associated with the DSMC grid. These 
particles are then implicitly mapped to the PIC nodes. For these 
reasons the PIC algorithm can be run in parallel by a simple tree sum 
of the charge density at the nodes. The overhead in this method is that 
each processor has a copy of the entire PIC grid and must calculate 
the electric field for the entire domain. However, these physical 
models only require a central differencing scheme for calculating 
the electric field. By making this simple change to the PIC part of 
the code, advantage is taken of the underlying parallelization of the 
DSMC code. For completeness, one time step of the algorithm is 
shown in Fig. 3. 

The behavior of the ions is of primary interest. A sufficient number 
of ion particles are needed in each cell to represent the ion density for 
the PIC method without introducing artificial gradients. Therefore, 
to increase the ion particle count without creating an overabun- 
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Fig. 3 One time step of the DSMC-P1C algorithm. 


dance of neutrals, a particle weighting scheme is used. The method 
outlined by Bird for species weighting only conserves momentum 
and energy on the average, not in each collision. 1 This approach 
is not conducive to flows with few'collisions. A method developed 
by Boyd 19 splits the particle with the larger weight into two parts. 
One of these has its momentum changed according to the collision 
dynamics, and the two parts are recombined afterwards. However, 
for collision pairs with vastly different velocities, this technique 
fails to represent the velocity distribution of the species with the 
larger weight. Charge-exchange reactions commonly involve col¬ 
lision pairs of this type. For this reason a new particle weighting 
scheme is used for the charge-exchange reactions, whereas Boyd’s 
scheme is used for momentum transfer. 

The present scheme is generalized for ion-neutral collisions re¬ 
gardless of which particle has the larger weight, but in these cases 
the neutral either has an equal or larger weight. The reaction can be 
described by the following equation: 

W„Xe, + W,XeJ -+ (W n - HMXe, + IV,Xe, + H^Xe* 

This scheme has the drawback that a third particle is created in each 
charge-exchange reaction. However, the mean free path is large 
enough (^1.2m at the thruster exit) that this is not a concern. Phys¬ 
ically, this shows that only W-, of the W n molecules undergo a col¬ 
lision and the others are unchanged. This leaves a distribution of 
W n — Wj molecules at the same velocity and the other W t with a 
much higher velocity. If these two particles were instead recom¬ 
bined using Boyd’s scheme, W n molecules would have a velocity 
between the extremes that none of these molecules physically would 
have. 

Most of the simulations employ 700,000-800,000 particles in 
1600 DSMC cells and 9500 PIC cells on an IBM SP2. These sim¬ 
ulations require about 50 CPU h distributed over four processors 
because the charge-exchange ion timescale is large compared to the 
propellent ions. Parallel efficiency ranged from 80 to 90%. The fluc¬ 
tuations are caused by the particle weighting scheme. Periodically, a 
load-balancing routine is used to distribute the particles more evenly 
among the processors. 

Flow Conditions 

The nominal operating conditions given in the experiments con¬ 
sidered are a current of 4.5 A, a discharge voltage of 300 V, a total 
flow rate of S.2-5.4 mg/s with 7% cathode split, and a facility back¬ 
pressure of 6 mPa. An ion temperature of 4 eV is assumed at the 
exit based on Ref. 11 for most of the cases and is compared with 
a 0.4 eV case (labeled case 3 in the figures). The electron temper¬ 
ature is assumed to be 3 eV based on Ref. 9 and is used to obtain 
the potential from the Boltzmann relation. These parameters define 
the input values, but it is necessary to separate them into densities 
and velocities for each species. The fraction of doubly charged ions 
is estimated as 25% by Manzella and Sankovic 10 for an SPT-100. 
This value along with the thruster’s efficiency and the assumptions 
about the neutrals allows exit plane properties for each species to 
be determined. 

The efficiency obtainable at nominal conditions is approximately 
0.5 for the SPT-100 Hall thruster 2 ^ A value of 0.54 is chosen 
for this investigation, which gives an ion* velocity close to King’s 
measurements. 12,13 The ions leave the thruster with a radial com¬ 
ponent to their velocity. 11 Oh and Hastings inverted an analytic 
function to choose a divergence angle for each particle. 7 The present 
study simplifies this process by assuming the divergence angle varies 
linearly with radial distance from the center of the annulus ring. A 
maximum of 10 deg is used in one case and 20 deg in another. A 
uniform ion velocity (0 m/s radial velocity) case is also considered 
for comparison with current density and heat flux measurements. 
Each generated ion is chosen from a Maxwellian distribution with 
the specified average velocity. These velocity profiles are combined 
with a uniform ion density. See Table 1 for a list of thmster exit 
conditions used for the sensitivity study. 

Table 2 shows the input conditions for examining the effects of 
a varying electron temperature. The main difference from the other 
input conditions is that these cases have only 10% of the propel¬ 
lent ions assumed to be doubly charged, 13 but this negligibly affects 
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Table 1 Ion conditions at thruster exit „ 
for sensitivity study 



Speed, 

Velocity 

Temperature, 

Case 

km/s 

variation 

eV 

1 

17.6 

/(r): 10 deg 

4 

2 

17.8 

/(r): 20 deg 

4 

3 

17.6 

/(r): 10 deg 

0.4 

4 

17.5 

Uniform 

4 

Table 2 

Input conditions for modeling study 

Case 

7V, eV 

T e profile 

n e oo, nr 3 

Ref* 

3.0 

Constant 

6.5*15 

VT1 

1.0 

Varying 

6.5*15 

VT2 

2.0 

Varying 

6.5*15 

VT 3* 

2.0 

Varying 

1.3*15 




Fig. 5 Comparisons of current density at radial distances of 50 cm and 
1 m. 



Fig. 6 Comparisons of current density at a radial distance of 50 cm. 


aggregate ion properties in the plume. The value for reference elec¬ 
tron number density determines the magnitude of the potential, but 
has little effect on the differences in the potential. The value of the 
electron temperature for the variable electron temperature cases rep¬ 
resents the far-field value imposed. The cases marked with asterisks 
include a full chamber geometry simulation as well as a smaller 
domain one. 

Results 

Measurements of radial velocity by Manzella 11 indicate an almost 
linear variation with radial position 11 mm from the thruster exit. A 
comparison of the simulations with this data is shown in Fig. 4. The 
simulations also give an almost linear variation that is reasonably 
close to the data at this axial location.* There are significant discrep¬ 
ancies between simulation and data. Neither case reproduces the 
asymmetry found experimentally because it was not included in the 
thruster exit plane profile for this investigation. 

Comparisons of current density for various simulations with ex¬ 
periments are presented at two different radial locations in Figs. 5-7. 
The base case (case 1) shown in Fig. 5 agrees well with King’s 
experiment. 12 The comparisons in Figs. 6 and 7 indicate theinsen- 
sitivity of the computed current density to the variations in inlet ion 
profiles. The variations are more pronounced close to the center¬ 
line; however, each case is within a factor of about two of the data. 
Comparing cases 1,2, and 4 indicates the significance of the diver¬ 
gence angle. Clearly, the higher angle leads to a lower peak value at 
these radial distances. Experiments show a more defined peak than 
the simulations. The low ion temperature case (case 3) gives a more 
defined peak, but an underprediction of current density occurs at 
higher angles. 

In Figs. 8-10 the current density is separated into density and 
velocity components. The velocity comparisons in Figs. 8 and 9 
indicate agreement within experimental uncertainty up to about 
50 deg. The simulation results only present the average ion velocity 



Fig. 7 Comparisons of current density at a radial distance of 1 m. 

of the propellent ions, neglecting the charge-exchange ions. The ex¬ 
perimental data were obtained from a retarding potential analyzer 
(RPA). 13 Numerical differentiation of the current vs voltage data 
approximates the derivative of this curve, which is proportional to 
the energy-per-charge distribution function. For singly charged ions 
of uniform mass, the energy distribution function also represents 
the velocity distribution function. By normalizing and numerically 
integrating this velocity distribution function, the average velocity 
is obtained. This process is likely to minimize the effects of the 
charge-exchange ions because they only contribute to the current 
measurements at low voltages. Also, the current scales with the flux 
of the charge-exchange ions. Thus, the density that is needed to 
obtain the actual contribution of the charge-exchange ions to the av¬ 
erage velocity may be obscured. Therefore, comparing the average 
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Fig. 9 Comparisons of average ion velocity at a radial distance of 1 m. 



Fig. 10 Comparisons of ion density at a radial distance of 50 cm. 


ion velocity of only the propellent ions from the simulation with 
the experimental measurements is reasonable. The agreement is en¬ 
couraging, The ion velocity and current density are used to obtain 
the experimental number density shown in Fig. 10. The simulations 
allow the ion density to be given directly. Both total density and the 
propellent ion density of the simulations are shown in Fig. 10. The 
propellent ion density agrees well at low angles, but it is too low 
at higher angles where the charge-exchange ions represent a higher 
fraction. The profile of the total density for case 1 compares favor¬ 
ably with the experiment, but the magnitude is too high. Following 
the same reasoning as described earlier, the measured ion density 
should be between the total and propellent densities from the simu¬ 
lation. In contrast to the ion velocity comparisons in Figs. 8 and 9, 
case 2 does not represent the ion density as well as case 1. 


Heat flux measurements were obtained by a heat flux probe 
consisting of two Schmidt-Boelter thermophile transducers. 13 The 
probe was biased to ground potential, which was between 5 and 8 V 
below the plasma potential, to eliminate electron current heating. 
Thus, the measured heat flux includes both neutrals and ions. The 
simulations give the contributions of the ions and the propellent 
neutrals. The background neutrals contribution is negligible. It is 
estimated to be on the order of 1 W/m 2 . In Figs. 11-13 compar¬ 
isons of the simulations with heat flux measurements are shown. 
There is agreement near the thruster centerline at 50 cm. At 1 m the 
simulated heat flux is within a factor of three of the measurements. 
At large angles from the centerline, where charge-exchange ions 
dominate, there is an underprediction by two orders of magnitude at 
both locations. The charge-exchange ions have much lower veloci¬ 
ties than the propellent ions, and the heat flux scales with velocity 



Fig. 11 Comparisons of heat flux at a radial distance of 50 cm. 



Fig. 12 Comparisons of heat flux at a radial distance of 1 m. 



Fig. 13 Comparisons of total heat flux at a radial distance of 50 cm. 
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Table 3 Integrated heat flux 
at various radial locations 


Case 

50 cm, W 

1 m, W 

King 

1692 

2339 

1 

630.0 

449.4 

2 

637.8 

460.0 

3 

609.6 

410.1 

4 

618.3 

424.8 

Ref 

652.5 

454.9 

VT 1 

678.0 

468.6 

VT2 

678.5 

472.6 

FG VT 

680.3 

478.3 



Fig. 14 Comparisons of electron density with measurements by 
Myers 9 at a radial distance of 31 cm. 

cubed. The spreading of the propellent ions does not improve the 
agreement at the higher angles significantly, and it leads to poorer 
agreement near the centerline (Figs. 11 and 12). As with the ion den¬ 
sity comparisons, a divergence angle ranging from —20 to 20 deg 
does not lead to as good agreement as the other cases. The ion tem¬ 
perature has a surprisingly small effect on the heat flux (Fig. 11). The 
ion temperature affects the spreading of the beam without altering 
the magnitude of the mean velocity. The higher temperature allows 
more spreading. The heat flux at higher angles should increase with 
ion temperature because its magnitude varies as velocity cubed. To 
estimate the effect on the ions of biasing the probe, an energy of 
8 eV is added per ion to the heat flux calculations and is labeled q m 
in Fig. 13. Clearly, it is only significant at the larger angles, but it 
does not make a substantial difference. Upon integrating the heat 
flux from the experiments, it is found that this heat transfer rate is 
greater than the power put into the system. Integrated values for the 
heat flux are shown in Table 3. The electrical power is only 1350 W. 
The simulations give values indicative of the measured efficiencies 
(?«50%) at a radial distance of 50 cm. By 1 m collisions with the 
background neutrals have decreased the energy of the ion flow. The 
simulations, which include a varying electron temperature model, 
are also included in the table, but plots of the heat flux show only 
moderate differences from the reference case. 

Calculated values for electron density are compared with mea¬ 
surements by Myers and Manzella 9 at a radial distance of 31 cm in 
Fig. 14. Cases 1 and 2 show the best agreement. The other cases 
have inlet conditions that lead to less spreading of the ions. The elec¬ 
tron densities from the other set of test cases are shown in Fig. 15. 
Including a varying electron temperature affects the spreading of 
the ions. The radial electric field is larger than it would be without a 
temperature gradient. The case labeled “Ref’ is nearly identical to 
case 1. It is clear from Fig. 15 that the case labeled “VT 3“ agrees 
best with the measurements. 

Comparisons with measurements of current density in the plume 
near field by Kim 14 are shown in Figs. 16 and 17 at three axial 
locations. At an axial distance of 10 mm, the agreement is reasonable 
regardless of varying the electron temperature. However, by 100 mm 
the effects of this temperature gradient are more apparent. The ion 



Fig. 15 Comparisons of electron density with measurements by 
Myers 9 at a radial distance of 31 cm. 




Fig. 17 Comparisons of current density from variable electron tem¬ 
perature case with measurements by Kim. 14 

beam reaches the axis more quickly than is predicted by the reference 
case. The variable electron temperature case more closely matches 
the shape of the measurements at 100 mm than the reference case. 
Its peak is about 50% higher than that of the reference case. Even by 
50 mm the current density is higher at the axis for a varying electron 
temperature. For these near-field comparisons there is negligible 
difference between the different variable electron temperature cases; 
therefore, only the first is shown. Integrated values are compared in 
Table 4. At 10 and 100 mm, varying the electron temperature gives 
only a 4% difference in integrated current. 
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Table 4 Integrated cu rrent at various axial locations 

Case 10 mm 25 mm 50 mm 100 min 

Ktei 3.97 A 4.92 A 4.95A 4.51 A 

Ref 4.04 A - 4.09 A 4.19 A 



Radial Distance (m) 


Fig. 18 Comparisons of ion density for the various simulations at an 
axial distance 4 cm behind the thruster. 



Fig. 19 Comparisons of axial current density for the various simula¬ 
tions at an axial distance 4 cm behind the thruster. 

For contamination concerns the backflow and far-field regions 
are more significant than the near field. The next few figures com¬ 
pare the reference case with the variable electron temperature cases. 
Figure 18 shows the densities of the various ions 4 cm behind the 
thruster. Clearly, the charge-exchange ions are most dominant. This 
is also true well above the thruster in the radial direction because the 
propellent ions are predominantly axially directed. Comparisons of 
the charge-exchange ions both behind the thruster and above it for 
the various cases indicate the significance of the magnitude of the 
far-field electron temperature. More importantly, the current density 
shows differences at these locations (Figs. 19 and 20). Behind the 
thruster the reference case has a magnitude that is a factor of two 
higher than the TV = 1 e V case near the thruster and 50% higher near 
the top of the domain. The results are similar at 70 cm above the 
thruster as well in both the axial and radial directions. 

Simulations, which include the full vacuum chamber geometry, 
allow comparison with experimental data at large angles in the back- 
flow region. The current density from the full geometry cases are 
compared with measurements by King 13 at a radius of 50 cm and 
1 m in Fig. 21. Even at large angles behind the thruster, good agree¬ 
ment is obtained for both models. The relative insensitivity of the 
current density to the varying electron temperature coupled with 
the differences in electron number density suggests that the far-field 
electron temperature is important. 



Fig. 20 Comparisons of current density for the various simulations at 
a radial distance 0.7 m above the thruster. 




Fig. 22 Comparisons of potential with measurements by Marrese 
(from private communication with A. Gallimore) at an axial distance of 
48 cm. 

The plasma potential is compared with measurements taken by 
Marrese and Gallimore (C. Marrese and A. Gallimore, University of 
Michigan, May 1999, private communication) at an axial distance 
of 48 cm in Fig. 22. In the simulations the potential drops below 
zero because of the choice of the reference electron density used in 
inverting the Boltzmann relation. The experimental measurements 
here and by King and Gallimore 13 (in Fig. 23) do not become neg¬ 
ative. The third variable electron temperature simulation is shown, 
which uses a lower value for the reference number density based on 
Marrese’s measurements. This value also agrees with the far-field 
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Fig. 23 Comparisons of plasma potential with measurements by 
King 13 at a radial distance of 50 cm. 


electron density of the simulations. This simulation matches Mar- 
rese’s data well. The first two variable electron temperature cases 
have the same shape as Marrese’s data but are negative like the ref¬ 
erence case. The agreement with King’s data (in Fig. 23) is less 
encouraging. The shape of the potential profile measured experi¬ 
mentally is explained in Ref. 21 through a detailed analysis of the 
effect on plasma potential of the magnetic configuration of several 
different Hall thrusters. 

Conclusions 

A computer code has been developed for simulating Hall thruster 
plumes. In this study the focus was on examining the sensitivity of 
the flowfield to various ion inlet conditions as well as to different 
models for computing the electric fields based on the electron mo¬ 
mentum equation. Simulations were performed for the SPT-100 at 
nominal operating conditions, and comparisons with measured data 
were presented for a number of plume properties including current 
density in the near-field, far-field, and backflow regions, and for 
electron density, heat flux, and plasma potential. 

Numerical issues such as computational time and memory cause 
restraints on the ability to perform large-scale simulations. A particle 
weighting scheme was used to allow a sufficient number of ions 
in the computational domain without an overabundance of neutral 
particles. This scheme better represents the velocity distribution than 
a weight-by-species method. The parallel nature of the code allows 
large-scale simulations. These two numerical techniques offer the 
capabilities of accurately simulating the physical models. 

Calculations of far-field ion properties are affected by variations 
in the ion inlet profile. However, the comparisons with experimental 
data show only moderate differences. Each of the various simula¬ 
tions agreed well with experimental measurements of far-field ion 
current density. 

Including a varying electron temperature was found to be signif¬ 
icant in the modeling of Hall thruster plumes. The electron temper¬ 
ature gradient significantly affects the near-field of the plume. Sim¬ 
ulations incorporating the varying electron temperature model gave 
better agreement with the near-field current density measurements 
than those using the isothermal Boltzmann relation. Results suggest 
that outside of this region the plasma behavior is not dictated by 
this near-field behavior, but rather by the value of the electron tem¬ 
perature. Thus, the far-field electron temperature is also important. 
This temperature is likely to be lower than the near-field value. The 
simulations agreed well with the measurements of electron num¬ 
ber density. Agreement with plasma potential was dependent on the 
value for the reference density used in the physical models. 

Simulating the full geometry of the experimental facility is impor¬ 
tant to consider measurements for code validation of the backflow 


region. This can also aid in determining any possible effects from the 
facility. The simulations were found to represent the current density 
very well. Even at large angles behind the thruster, good agreement 
was obtained with experiment. 
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